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The nature of blood as a suspension of red blood cells makes computational hemo- 
dynamics a demanding task. Our coarse-grained blood model, which builds on a 
lattice Boltzmann method for soft particle suspensions, enables the study of the 
collective behavior of the order of 10 6 cells in suspension. After demonstrating the 
viscosity measurement in Kolmogorov flow, we focus on the statistical analysis of 
the cell orientation and rotation in Couette flow. We quantify the average incli- 
nation with respect to the flow and the nematic order as a function of shear rate 
and hematocrit. We further record the distribution of rotation periods around the 
vorticity direction and find a pronounced peak in the vicinity of the theoretical 
value for free model cells even though cell-cell interactions manifest themselves in 
a substantial width of the distribution. 
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1. Introduction 

Human blood can be approximated as a suspension of deformable red blood cells 
(RBCs, also called erythrocytes) in a Newtonian liquid, the blood plasma. The 
other constituents like leukocytes and thrombocytes can be neglected due to their 
low volume concentrations. Typical physiological RBC concentrations are around 
50 %. In the absence of external stresses,erythrocytes assume the shape of biconcave 
discs of approximately 8 (jjxi diameter [l( ■ An understanding of their effect on the 
rheology and the clotting behavior of blood is necessary for the study of pathological 
deviations in the body and the design of microfluidic devices for improved blood 
analysis. 

Well-established methods for the computer simulation of blood flow either con- 
sist of an elaborate model of deformable cells 0, Q or restrict themselves to a con- 
tinuous description at larger scales Q. Our motivation is to bridge the gap between 
both classes of models by an intermediate approach: we keep the particulate nature 



Article submitted to Royal Society 



TfcjX Paper 



2 



F. Janoschek, F. Mancini, J. Harting, and F. Toschi 



of blood but simplify the description of each cell as far as possible. The main idea 
of the model [B| is to distinguish between the long-range hydrodynamic coupling 
of cells and their complex short-range interactions. Our method is well suited for 
the implementation of complex boundary conditions and an efficient parallelization 
on parallel supercomputers. This is important for the accumulation of statistically 
relevant data to link bulk properties, for example the effective viscosity, to phe- 
nomena at the level of single erythrocytes. In this article, we present an alternative 
method of viscosity measurement and a study of the orientational dynamics of cells 
in a sheared suspensions. The improved understanding of the dynamic behavior of 
blood might be used for the optimization of macroscopic simulation methods. 

2. Coarse-grained model for blood flow simulations 

We apply a D3Q19/BGK lattice Boltzmann method for modeling the blood plasma 
The single particle distribution function n r (x, t) resembles the fluid traveling 
with one of r = 1, . . . , 19 discrete velocities c r at the three-dimensional lattice 
position x and discrete time t. Its evolution in time is determined by the lattice 
Boltzmann equation 

n r (x, t) — n° q (p(x, t), u(x, t)) , 
n r {x + c r> t + l) = n r (x,t)-Cl with Q= v ' ' r VFV ' ; ' v 1 " (2.1) 

T 

being the BGK-collision term with a single relaxation time r. The equilibrium dis- 
tribution function n° q (p, u) is an expansion of the Maxwell-Boltzmann distribution. 
p(x, t) = J2 r n r{ x i t) and p(x, t)u(x, t) = ^ r ?v(x, t)c r can be identified as density 
and momentum. In the limit of small velocities and lattice spacings the Navier- 
Stokes equations are recovered with a kinematic viscosity of v = (2r — l)/6, where 
r = 1 in this study. 

For a coarse-grained description of the hydrodynamic interaction of cells and 
blood plasma, a method similar to the one by Aidun et al. modeling rigid particles of 
finite size is applied Q- It can be seen as a bounce-back boundary condition which 
is enhanced by a correction term C = 2a Cr p(x. + c r , t) c r v/cg that accounts for a 
possible local boundary velocity v. The lattice weights a Cr and the speed of sound 
c s are constants for a given set of discrete velocities. The resulting bounce-back rule 

n r (x + c r ,t + 1) = ni(x + c r ,t) +C with n*(x, t) = n r (x, t) - fi (2.2) 

and f specifying the direction opposite to r is applied to all fluid distributions 
that according to Eq. 12.11 would travel along a link that crosses the theoretical cell 
surface. The momentum Apf p = {2rif + C) c f which is transferred during each time 
step by each single bounce-back process is used to calculate the resulting force on 
the boundary. 

Instead of the biconcave equilibrium shape of physiological RBCs we choose a 
simplified ellipsoidal geometry that is defined by two distinct half-axes R\\ and R± 
parallel and perpendicular to the unit vector 6 2 ; which points along the direction 
of the axis of rotational symmetry of each cell i. Since the resulting volumes are 
rigid we allow them to overlap in order to account for the deformability of real 
erythrocytes. We assume a pair of mutual forces F pp = ±2n° q (p, u = 0)c r at 
each cell-cell link. This is exactly the momentum transfer during one time step 
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Figure 1. Outline of the model. Shown are two cells with their axes of rotational symmetry 
depicted by vectors. The volumes defined by the cell-cell interaction are approximately 

ellipsoidal ( — ). The smaller ellipsoidal volumes ( ) of the cell-plasma interaction are 

discretized on the underlying lattice. (Online version in color.) 

due to the rigid-particle algorithm for a resting particle and an adjacent site with 
resting fluid at equilibrium and initial density p. While these surface forces solely 
serve to compensate the lack of fluid pressure inside the particles, we additionally 
implement phenomenological pair potentials for neighboring cells that model the 
volume exclusion of soft and anisotropic particles in an effective way. For the sake 
of simplicity, we use the repulsive branch of a Hookian spring potential 
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for the scalar displacement r.^ of two cells i and j. With respect to the disc-shape 
of RBCs, we follow the approach of Berne and Pechukas Q and choose the energy 
and range parameters 



e(6i,6j) = e 1 - x 2 {biOjY 



-1/2 
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as functions of the orientations 6^ and bj of the cells and their normalized center 
displacement f^. We achieve an anisotropic potential with a zero-energy surface 
that is approximately that of ellipsoidal discs. Their half-axes parallel cry and per- 
pendicular £7j_ to the symmetry axis enter Eq. 12.41 and Eq. 12.51 via a = 2cr± and 
X = (o" 2 — cr 2 ^)/(c| + <j\) whereas e determines the potential strength. 

Fig. [T] shows an outline of the model. Two cells i and j surrounded by blood 
plasma are displayed. For clarity, the three-dimensional model is visualized in a two- 
dimensional cut. Depicted are the cell shapes defined by the zero-energy surface of 
the cell-cell potential Eq. 12.31 with Eq. 12.41 and Eq. 12.51 that can be approximated 
by ellipsoids with the size parameters fry and a± as half axes. The cells are free to 
assume continuous positions and orientations and Oj. In consequence, also the 
center displacement vectors and r^ x between cells are continuous. Still, for the 
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cell-plasma interaction an ellipsoidal volume with half axes Ru and R±_ is discretized 
on the underlying lattice. For additional clarity, the discretization is drawn only for 
cell j. The forces and torques emerging from the interaction of the cells with other 
RBCs and the fluid are integrated by a classical molecular dynamics code in order to 
evolve the system in time. The conversion from lattice units to physical units is done 
by multiplication with 8x = 2/3 ^m, St = 6.80 x 1CT 8 s, and 8m = 3.05 x 10~ 16 kg 
for length, time, and mass respectively. As a convention, primed variables are used 
whenever we refer to quantities specified in physical units. The model parameters 
are chosen as R' ± = 11/3 /im, R'« = 11/9 fim, a' ± = 4 /im, ai = 4/3 ^m, and 
e' = 1.47 x 10~ 15 J. For more detailed information see 



3. Results 

(a) Apparent bulk viscosity measurement in Kolmogorov flow 

In our previous work, the apparent viscosity of the bulk of the suspension was 
measured in simulations of unbounded Couette flow [B| . Though this setting is easy 
to comprise analytically, its efficient and precise implementation proves a non-trivial 
challenge since the suspended particles would have to be enabled to stretch across 
the sheared boundary. We therefore demonstrate an alternative method to deter- 
mine the apparent viscosity that is compatible with completely periodic boundaries 
and is based on so-called Kolmogorov flow: a sinusoidally modulated shear flow. We 
proceed as in Benzi et al. Q and apply to the full suspension a sinusoidal body 
force 

f z (x) =/osin[fc(z-l/2)] (3.1) 

pointing into the z-direction. /o is the amplitude and k the wave number in x- 
direction. At steady state, the spatial variation of the shear stress 

d x o- xz {x) = f z (x) (3.2) 

matches the external forcing. Integration of Eg. 13.21 together with Eg. 13.11 yields an 
analytic expression for the shear stress 

g xz {x) = cos [k{x - 1/2)] . (3.3) 

The y- and z-averaged local shear rate can be evaluated numerically as 

ixz{x) = {d x u z {x)) y z . (3.4) 

After equilibration, a simulation of N x x N y x N z lattice sites results in N x numbers 
for the apparent viscosity 

/ \ °~ X z{x) /„ rl 

Ma PP (a;) = - — j-t (3.5) 

^ X z \X) 

covering a range of shear rates that depends on /o and k. 

Compared to a viscosity measurement in Couette flow the above proce- 
dure Q has the benefit of taking advantage of the whole simulation volume since 
there are no possibly unphysical boundary regions. Additionally, each single mea- 
surement yields data for many different shear rates. On the other side, the non- 
constant shear stress causes inhomogeneities in the local cell volume concentration 
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Figure 2. Comparison of viscosity data retrieved from Kolmogorov flow at different wave 
numbers k and amplitudes of forcing /□ after 68.0 ms of equilibration with data measured 
in Couette flow as in The hematocrit is $* =43% with A$ = 0.5% while the error 
bars were drawn according to the standard deviation after binning and averaging of the 
original data. (Online version in color.) 

which are much more pronounced than in the case of Couette flow. Therefore, 
we disregard all viscosity data (x,^ app (x)) for which $(x) [$* - A$, $* + A$] 
with $* being the volume concentration of interest and 2A$ a small interval of 
tolerance. Fig. [2] compares A'appd'Txzl) resulting from such measurements at vary- 
ing /o and N x with k = 2n/N x after 68.0 ms with Couette flow measurements as 
described in @ (N x = 128). In the plot, $* = 43% and A$ = 0.5%. We further 
average the data within bins with a width of A[ln(7^ s)] = 0.5 in order to reduce 
statistical noise. The good agreement proves the feasibility of the method. 

(6) Statistical analysis of cell orientations in Couette flow 

Due to its applicability to huge numbers of particles and the fact that cell orien- 
tations are accessible directly in the model, our simulation method is particularly 
suited for the statistical analysis of RBC orientations at high volume concentra- 
tions. In the following, we simulate Couette flow as described in However, with 
a size of N x — N z =512 and N y = 64 lattice sites, the volume is considerably larger 
and contains about 3 x 10 4 cells at a hematocrit of $ = 45 %. The velocity gradient 
is aligned with x- and the flow with the z-axis. To exclude boundary effects in the x 
direction, we restrict ourselves to the analysis of those cells which stay at a lateral 
position 512/3 < x < 2 x 512/3 during the whole simulation. The visualization 
of the cell volumes defined by the model potential is shown for a smaller system 
in the inset of Fig. [3] It is clear that also at physiological volume concentrations 
there is a preferential orientation of cells. In order to quantify this observation, the 
orientation vector of every cell i is multiplied with —1 where necessary to achieve 
(6i) x > 0. The angle 9 is measured between the positive x-axis and the resulting 
normalized orientation vector. Fig. [3] displays the most probable value 6* that is 
determined by fitting a Gaussian to the distribution of angles. For $ = 45 %, com- 
pared to 9 — which would mean alignment parallel to the flow, there always is a 
positive inclination 9* which is directed against the vorticity of the shear (see inset 
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Figure 3. Most probable value of the angle 9 between the axis of symmetry of a model cell 
and the rr-direction as a function of the shear rate |7i 2 | at <J> = 45 %. The statistical error 
of 9 is smaller than the symbols. The inset shows alignment of cells due to shear flow in 
a smaller system. (Online version in color.) 

of Fig. [3]). When increasing the shear rate from around 300s -1 to 6000 s -1 , this 
inclination becomes considerably smaller: 8* decreases from 14.5° to 3.4° . 

Studying 8 allows to quantify the orientation with respect to the direction of 
the flow. To quantify the actual amount of orientational order, the nematic order 
parameter known from liquid crystal physics (see for example (lo| ) is a better 
choice. This parameter is usually defined as the largest eigenvalue A + of the nematic 
order tensor Ski = (3ofe6; — Ski),: t /2 which is obtained as the average over different 
time steps t and cells i within the volume of interest. Possible values for A+ are 
comprised between and 1 indicating complete disorder and order. Fig. U depicts 
A+ for different shear rates. A decrease of nematic order with increasing shear rate 
is visible. The inset shows the dependency on the hematocrit $ at a fixed shear 
rate of y xz = (3.0 ± 0.1) x 10 3 s _1 . For comparison, also A+ resulting for a single 
ellipsoidal particle that tumbles with 6 perpendicular to the vorticity direction is 
shown. Apparently, the higher nematic order at $ = 45 % is caused by hindered 
tumbling motion while the reason for the lower A+ at $ = 11% is the lack of 
alignment in the vorticity plane. 

(c) Statistical analysis of cell rotations in Couette flow 

The above conclusions make clear that preferential orientations result from the 
averaged tumbling of many cells. The time evolution of the continuous tumbling 
angle 8 is plotted for an arbitrary selection of cells in the inset of Fig. [5j The respec- 
tive volume concentration is $ = 45 % and the shear rate j' xz = 6.1 x 10 3 s _1 . The 
cells keep a largely constant alignment for varying periods of time and, occasionally, 
flip over by an angle of tt along the vorticity direction. Thus, the angular velocity is 
strongly time-dependent. The probability distribution function of the time required 
for a rotation by tt which is measured as the time T between two flipping events is 
shown in the main part of Fig. [3] For decreasing volume concentrations, the distri- 
bution becomes narrower and its peak is shifted towards shorter times. However, 
with a width comparable to its average value, even at <J> = 11 % the distribution 
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Figure 5. Probability density functions of the time T' between two consecutive events 
of flipping by an angle of tt around the vorticity direction at different values of $ and 
i'xz ~ (3-0 ± 0.1) x 10 3 s~ 1 . As a reference, the vertical lines indicate the time required 
for a rotation around n according to the analytical solution for a single ellipsoidal particle 
at the respective shear rates [TJ. The inset visualizes the continuous tumbling angle 6 of 
selected cells as a function of time t' at $ = 45 %. (Online version in color.) 



deviates substantially from a delta peak. Therefore, even though typical tumbling 
periods in suspension do not differ much from the value for a freely tumbling ellip- 



soid [ll|, the effect of cell-cell interactions on our observables is significant for all 
<f> studied. 



4. Conclusion 

Our mesoscale model for blood Q was applied to bulk flow of up to 3 x 10 4 cells. 
We demonstrated the feasibility of viscosity measurement in Kolmogorov flow and 
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studied the statistics of the rotational cell positions and velocities at different values 
of the hematocrit. These observables are extendable also to more elaborate (e. g. 
deformable) cell models. Their quantitative study and measurement is a key to 
build more reliable cell based or continuum descriptions for blood. 
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